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Abstract 

Elastic moduli and dislocation core energy of the triangular solid of hard 
disks of diameter a are obtained in the limit of vanishing dislocation- an- 
tidislocation pair density, from Monte Carlo simulations which incorporates 
a constraint, namely that all moves altering the local connectivity away from 
that of the ideal triangular lattice are rejected. In this limit, we show that 
the solid is stable against all other fluctuations at least upto densities as low 
as /3cr^ = 0.88. Our system does not show any phase transition so diverging 
correlation lengths leading to finite size effects and slow relaxations do not 
exist. The dislocation pair formation probability is estimated from the frac- 
tion of moves rejected due to the constraint which yields, in turn, the core 
energy Ec and the (bare) dislocation fugacity y. Using these quantities, we 
check the relative validity of first order and Kosterlitz -Thouless -Halperin - 
Nelson -Young (KTHNY) melting scenarios and obtain numerical estimates 
of the typical expected transition densities and pressures. We conclude that 
a KTHNY transition from the solid to a hexatic phase preempts the solid to 
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liquid first order transition in this system albeit by a very small margin, easily 
masked by crossover effects in unconstrained "brute - force" simulations with 
small number of particles. 

I. INTRODUCTION 

One of the first continuous systems to be studied by computer simulations [|l],H is the 
system of hard disks interacting with the two body potential, 

V{r) = oo,r < a (1) 
= 0, r > cr 

where, a (taken to be 1 in the rest of the paper) the hard disk diameter, sets the length scale 
for the system and the energy scale is set by fc^T = 1. Despite its simplicity 0], this system 
was shown to undergo a phase transition from solid to liquid as the density p was decreased. 
The nature of this phase transition, however, is still being debated. Early simulations 
always found strong first order transitions. As computational power increased the observed 
strength of the first order transition progressively decreased! Using sophisticated techniques 
Lee and Strandburg |^ and ZoUweg and Chester found evidence for, at best, a weak first 
order transition. A first order transition has also been predicted by theoretical approaches 
based on density functional theory [^]. On the other hand, recent simulations of hard 
disks by Jaster, using as many as = 65536 particles find evidence for a continuous, 
Kosterlitz -Thouless -Halperin -Nelson -Young (KTHNY) transition from liquid to a 
hexatic phase, with orientational but no translational order, at p = 0.899. Nothing could be 
ascertained, however, about the expected hexatic to the crystalline solid transition at higher 
densities because the computations became prohibitively expensive. The solid to hexatic 
melting transition was estimated to occur at a density pc > .91. A priori, it is difficult 
to assess why various simulations give contradicting results concerning the order of the 
transition. In this paper we take an approach, complementary to Jaster's, and investigate 
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the melting transition of the sohd phase. We show that the hard disk sohd is unstable 
to perturbations which attempt to produce free dislocations leading to a solid — > hexatic 
transition in accordance with KTHNY theory 0. Though this has been attempted in the 



past [|TO],|rT| , numerical difficulties, especially with regard to equilibration of defect degrees of 
freedom, makes this task highly challenging. We also show that this transition lies close to 
a, first order, solid to liquid melting line. We calculate quantitatively the relative positions 
of the first order and the KTHNY transitions in the parameter space for this system and 
explain why earlier simulations failed to arrive at a consensus. 

The coarse grained density of a crystalline solid can be expanded as, p(r) = J2g Pgc**^ "", 
where (G) is a reciprocal lattice vector. The order parameters pc are complex, pc = 
Ip^jgiu G^ and the displacement vector u is the deviation of an atom from the nearest perfect 
lattice point R. If fiuctuations of the amplitude of pc can be neglected then a solid can 
be described in terms of u alone — - the fundamental assumption of elasticity theory. The 
elastic Hamiltonian for hard disks is given by, 

F = -Pe+ + B/2el + (p + P)(e-/2 + 2e,,), (2) 

where B is the bulk modulus. The quantity pe// = p + -P is the "effective" shear modulus 
(the slope of the shear stress vs shear strain curve) and P is the pressure. The hard disk 
solid, being a purely repulsive system, is always under a uniform hydrostatic pressure -P(p) 
at any density p. The Lagrangian elastic strains are defined as, 

_1 / duj duj duj duk \ , s 

~ 2 \dR, OR, dRk OR, J ' ^ ^ 

where the indices i,j go over x and y and finally, e+ = exx + ^yy, and e_ = e^x — ^yy 

In general a solid possesses two types of excitations, "smooth" phonons and "singular" 
dislocations, respectively. Long wavelength phonons inhibit long range order in 2-d solids 
so that the intensity of a Bragg refiection peak Iq ~ e"^^"^ , where the Debye Waller factor 
Wg ~ G'^a?'^^ I (d — 2) [a is the lattice parameter and d the number of spatial dimensions) 
diverges and order parameter correlations decay algebraically — - an example of Quasi Long 



Ranged Order (QLRO). We know that singular excitations, like dislocations, can drive a 
QLRO disorder transition (where correlations decay exponentially). This situation has 
been analysed by the KTHNY theory 0. 

The KTHNY- theory is presented usually for a 2-d triangular solid under zero external 
stress. It is shown that the dimensionless Young's modulus of a two-dimensional solid, 

K = ^ ^ 

v/3p{l + /i/(A + /i)}' 

where /i and A are the Lame constants, depends on the fugacity of dislocation pairs, 
y = exp{—Ec), where is the core energy of the dislocation, and the "coarse -graining" 
length scale I. This dependence is expressed in the form of the following coupled differential 
equations (the recursion relations) for the renormalization of K and y: 

^ = 3vr,^es.[-/o(-) --/,(-)], (4) 

dy , K . ^ K ^ , K . 

where Jq and Ji are Bessel functions. The thermodynamic value is recovered by taking the 
limit / oo. 

We see in Fig. (1) that the trajectories in y-K plane can be classified in two classes, 
namely those for which ?/ ^ as / ^ oo (ordered phase) and those y oo I ^ oo 
(disordered phase). These two classes of flows are separated by lines called the separatrix. 
The transition temperature T^. (or p^) is given by the intersection of the separatrix with 
the line of initial conditions K{p,T) and y = exp{—Ec{K)) where E^ ~ cK/IGtt. At the 
transition point the flow follows the separatrix so that the renormalized K jumps from IGvr 
to at the transition. The ordered phase corresponds to the solid (no free dislocations) 
and the disordered phase is a phase where free dislocations proliferate. Proliferation of 
dislocations however does not produce a liquid, rather a liquid crystalline phase called a 
"hexatic" with quasi- long ranged (QLR) orientational order but short ranged positional 
order. A second K-T transition destroys QLR orientational order and takes the hexatic to 
the liquid phase by the proliferation of "disclinations" (scalar charges). Apart from Tc there 



are several universal predictions from KTHNY- theory, for example, the order parameter 
correlation length and susceptibility has essential singularities (~ e''* " ,t = T/T^ — 1) near 
Tc. All these predictions can, in principle, be checked in simulations [Q. 

Note that, in order to use the KTHNY- theory to study the solid- hexatic transition in 
hard disks we have to bear in mind that for the hard disk solid, which is always under a 
uniform hydrostatic pressure -P(p), the effective shear modulus /ie// has to be used in the 
definition of K. 

The KTHNY- theory predicts when a 2-D solid becomes unstable to the proliferation 
of dislocations. However, there is a second possibility. The free energy of the liquid may 
become higher than that of the stable solid at a density smaller than that where the hexatic 
phase is recorded. This leads to a first order transition and a jump in density at the liquid - 
solid coexistence pressure (for simulations in the NVT- ensemble) instead of an intermediate 
hexatic phase. Often it is very difficult to distinguish the two possibilities as the history 
of simulation studies of hard disks shows. This is further complicated by the fact that 
KTHNY theory also predicts that the specific heat, or equivalently, in the case of the hard 
disk system, the compressibility, shows a smooth bump leading to a near fiat region in the 
pressure -density diagram. In Fig. (2a) we show the conventional situation where the dotted 
line designates the often observed first order transition. In Fig. (2b) we show Jaster's results 
where it is seen that instead of a fiat region in the P-p- curve or a Maxwell loop usually 
associated with a first order transition one gets instead a smooth bending over to a state 
with a high compressibility. Finite size effects which would be present in the first- order case 
are negligible. This would indicate the presence of a liquid- hexatic transition. The question 
of solid to hexatic transition is still open. It is worth noting that detailed finite size scaling 
of orientational order in this system is not necessarily in contradiction to this result. 



Why do simulations of hard disk solids find it so difficult to see a solid- hexatic transition? 
One of the reasons is, of course, the divergence of the correlation length as the system 
approaches the transition so that one requires large systems. This is complicated by the fact 
that in order to obtain equilibrated values of the dislocation density (oc y) one also needs 



very large simulation times because in a high density solid the diffusion of defects is very 



slow . To illustrate this point we have attempted to calculate the defect density of a hard 
disk solid in a Monte Carlo simulation. We perform conventional Monte Carlo simulations 
in the NVT ensemble with an usual Metropolis updating scheme for = 3120 particles. We 
choose a single density p = .92; a sequence of initial states are then constructed by adding 
extra complete rows of atoms (thereby increasing the density to pi > .92) and removing an 
equal number of atoms from the bulk at random. In equilibrium, these extra vacancies in 
the bulk should diffuse out and the lattice parameter adjust to fill in the gap. After about 
one million Monte Carlo steps we calculate the number of five coordinated (77,5) and seven 
coordinated (^7) atoms. Since our system cannot have free vacancies (due to our choice of 
ensemble) we expect in equilibrium ns = nj. The simulations at each pi is repeated for ten 
realizations of the initial state. Our results are shown in Fig. (3). We see that 7^ 77,7 
(except for the trivial case of pi = .92), the difference growing with pi as expected, and 
the statistical errors are very large. We therefore conclude that even for a relatively small 
system of 3120 particles, the equilibration of defects (vacancies in this case) is an extremely 
slow process. So it should not come as a surprise that brute force simulations of the hard 
disk solid fail to produce the true equilibrium phase. 

It may also happen, on the other hand, that KTHNY- theory fails due to the following 
reasons. Firstly, elastic theory itself may fail near the transition, so that amplitude or long 
wave length phonon fluctuations may destabilise the solid producing a continuous transi- 
tion. Though remote, this possibility has nevertheless been discussed in the literature ||15|| . 
Secondly, perturbation theory in y may break down because is too small (i.e. y too large) 
at the transition. Saito [|16| and Strandburg showed using lattice discretized versions 
of a dislocation Hamiltonian that KTHNY perturbation theory breaks down if < 2.7 at 
the transition. In our simulations of the hard disk system we check both these possibilities 
as well as the possibility of a first order transition. 

In the next section we discuss our simulations together with our method for computing 
elastic constants and core energies. We use these inputs to check for a first order transition 



and a KTHNY- scenario in Section III. We summarise and conclude this work in Section IV. 



II. ELASTIC CONSTANTS AND DISLOCATION CORE ENERGIES FROM 

CONSTRAINED SIMULATIONS 

One way to circumvent the problem of large finite size effects and slow relaxation due 
to diverging correlation lengths is to simulate a system which is constrained to remain 
defect (dislocation) free and, as it turns out, without a phase transition. Relatively small 
systems simulated for short times therefore yields thermodynamically accurate data in this 
limit. Surprisingly, we show that using this data it is possible to predict the expected 
equilibrium behaviour of the unconstrained system. It is worth mentioning that with an 
approach similar in spirit to the one followed here, we have obtained excellent results for 



the Kosterlitz -Thouless transition in the two -dimensional planar rotor model [jl8|, which 
has served as an important model in the development of the KTHNY theory after the 
proofs of the low temperature susceptibility divergence in this model and the existence 



of phase transitions without local order parameters in general were given. 

We simulate = 3120 hard disks in an (almost) square box. We have also simulated two 
additional systems of = 2016 and = 4012 particles in order to look for residual finite 
size effects. Our algorithm follows closely the usual Metropolis scheme for simulating hard 
disks. The simulation is always started from a perfect triangular lattice which fits into our 
box - the size of the box determining the density. Once a regular MC move is about to be 
accepted, we perform a local Delaunay triangulation involving the moved disk and its nearest 
and next nearest neighbors. We compare the connectivity of this Delaunay triangulation 
with that of the reference lattice (a copy of the initial state) around the same particle. If any 
old bond is broken and a new bond formed (Fig. (4)) we reject the move since one can show 
that this is equivalent to a dislocation - antidislocation pair separated by one lattice constant 
involving dislocations of the smallest Burger's vector. Note that, (i) only dislocation pairs 
of smallest Burger's vector are eliminated, dislocations of higher topological charge cost 
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higher energy and may not be relevant at the densities where a melting transition is usually 
observed; (ii) other fluctuations e.g. long wavelength phonon fluctuations and fluctuations 
of the amplitude of the order parameter (spontaneous production of voids in the system) 
are not eliminated as long as they preserve connectivity. The fraction of moves p which 
are rejected because they violate the constraint is stored. Next, we need a method to 
calculate elastic constants accurately in our simulations making sure that we extrapolate to 
the thermodynamic limit. Such a method has been recently developed by us and discussed 
in detail elsewhere [^. Below we include a brief description for completeness. 



Since we have a dislocation free system, we can always associate an ideal, static, "ref- 
erence" lattice point R with every hard disk all through the simulation and calculate 
UR(t) = R'(i) — R- Microscopic strains ejj(R) can be calculated now for every reference 
lattice point R. Next, we coarse grain (average) the microscopic strains within a sub-box of 
size Lb, 



and calculate the {Lb dependent) quantities, 

S^X = < e+e+ > (5) 



'^33' — 4 < €xy€xy > 



The sub-blocks may be constructed by simply dividing the entire box of size L into integral 
number of smaller boxes, as done in this calculation so that L/L^ = an integer, or multiple 
sub-boxes of arbitrary size Lb < L can be constructed within the simulation cell, as in Ref. 
2T| . Lastly, quantities in the thermodynamic limit are obtained by fltting data to the form, 



C-tJfc _ coo 



where the index 7 = -|-, — , 3 the function, \1/(q;) is deflned as, 

2 . /•! 



^(a) = -o? [ [ dxdy Ko{aJx^ + y"^) 

TT Jo Jo ^ 



Kq is a Bessel function and ^ is the correlation length for the ee- correlations. 

The elastic constants in the thermodynamic limit are obtained from, the set: B = l/25'^_,_ 
and fieff = \/2S°^_ = 1/2S^. The last two equations for fiejf serve as a stringent internal 
consistency check and yields an accurate error estimate for this quantity. There are two 
ways to obtain the fluctuations S!^J^ for every sub-block size Lf, in Eq. (5). One can either 
accumulate < e^e^ > directly or construct histograms of the block strains and obtain 
Sjj by fitting Gaussian profiles to the normalized probability distributions of e-y for every 
block size Lb. Again this constitutes another excellent consistency check and a measure of 
the statistical uncertainties involved. We accumulate data till all these uncertainties are less 
than a percent. Residual finite size effects obtained by repeating the entire procedure for 
= 2016 and 4012 particles for a few densities are also seen to be within the same limit of 
accuracy. 

There are several distinct advantages of our method: In general our method works for any 
system for which instantaneous configurations can be obtained (for example either from other 
simulations or from real experiments). We obtain directly the finite size scaled results from 
a single simulation. As discussed above there are a number of stringent internal consistency 
checks which can be used to obtain very accurate data. In spite of this our method is easy to 
use and the computational complexity is not more than calculating for eg. pair correlation 
functions. This method can be easily adapted for calculating local strains and stresses in 
inhomogeneous situations. 



III. RESULTS AND DISCUSSION 

Our results for the elastic moduli, the pressure and the fraction of moves, p, rejected due 
to the topological constraint discussed above are given in Table I as a function of density. In 
Fig. (5), we compare our results for the bulk and shear moduli with the data of two previous 
simulations of Ref . |T^ and Ref . ||2^ . We also compare our simulation results to estimates 



from free volume theory p2[ in the simplest, independent cell approximation. Within this 
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approach the Helmholtz free energy per particle is given by / = log(f/), where the available 
free volume, vj = (a — 1)^/ pc and the close packed density pc = 2/ v^. Other thermodynamic 
quantities can be obtained by successive differentiation, viz. 

P = (6) 

X — 1 

B = P{1 



2{x - 1)' 
B 

Peff - ^ 

where x = \J~Pcfp and we have used the Cauchy relation, strictly valid only for a harmonic 
solid for our estimate of the effective shear modulus Peff- Note that the free volume 



elastic moduli and the pressure diverge |2j as p ^ pc- 

We see that our bulk modulus interpolates smoothly from the free volume values at high 



densities to those of Ref. at low densities. Overall, the differences between the three 
sets of data are small. Our values for the shear modulus agrees well with the free volume 
results at high density, but at low densities they are smaller than all other estimates though 
close to those of Ref. |T0 |. 



Once the elastic constants are obtained we can analyze in detail the two competing 
scenarios viz. first order solid- liquid transition or KTHNY- transition to the hexatic phase. 



A. Equation of state, free energy and first order melting 

First of all, we should point out that our constrained simulations allow us to obtain 
elastic constants up to a density as low as p = .88, far below the density p = .899 where the 
transition to the liquid is expected to occur 0|, which implies that amplitude and phonon 
fluctuations cannot destabilize the solid. So an ordinary second order transition is ruled out. 
However, there can always be a first order transition if the free energy of the liquid becomes 
lower than that of the perfect solid. 

In order to investigate this question we obtain the equation of state -P(p) and the Gibbs 
free energy g{P) of the liquid and the solid. 
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To obtain the equation of state of the hquid we use the semi -empirical, accurate, ana- 
lytical form by Santos et. al. ||2^, which is in excellent agreement with computer simulation 
data . The pressure is given by, 

P/p = {l-2r^ + ^^ryT^ (7) 

Tc 

where the packing fraction 77 = (7r/4)p and rjc is the packing fraction at close packing. The 
Helmholtz free energy per particle, 

/(p)= Pdr^'^l^ + U (8) 
Jo r]' 

where the ideal gas Helmholtz free energy per particle fid = log(p) — 1. The Gibbs free 
energy g{P) is then obtained by the standard Legendre transformation, g = f + P/ p. In 
addition we use the data of Jaster ^ in the transition region to obtain a revised estimate of 
the free energy. This is done by fitting Jaster's data for P{p) to a polynomial for p > 0.85 



which matches the results of Santos et. al. ||23[ for p < 0.85. From this equation of state we 



can obtain the Helmholtz and hence the Gibbs free energy by integrating starting from the 
value given by Eq.(7) at p = 0.85. 

The equation of state for the solid is obtained by integrating our bulk modulus values 
using the result of Bladon and Frenkel at p = 1.049 as the reference pressure (P = 22.00). 
The Gibbs free energy is obtained by further integration again using the result obtained for 
the free energy in Ref. |2^ at p = 1.049 as a reference {g = 25.64). 

The possible (first order) transitions can be located by equating the Gibbs free energies. 
The slope discontinuity gives the (inverse) density difference of coexisting phases. We find 
immediately, that all the free energies have very similar slopes (see Fig. 6) so that any 
possible first order transition would have only a small jump in the density. It also implies 
that small errors in the free energy of our reference state makes a large difference in the co- 
existence pressure. We have therefore reduced the reference free energy by a small amount 
(< 4%) so that the coexisting pressure Pi = 9.2 — the value found in most recent simulations 
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Using the semi-empirical free energy of Santos et. al. we obtain a (meta stable) 



first order transition with pi = 0.871 and = 0.912 as observed in early simulations 
Of course, this estimate of pi is only a lower bound, as the theory of Ref. is expected 
to overestimate the free energy. The free energy from Jaster's data is lower and almost 
completely parallel to that of the solid suggesting a very weak first order transition if at all. 
In this case we get a slope difference < 1.3% (viz. pi = 0.899 and ps = 0.911) - well within 
our numerical accuracy (Fig. 6). 

B. Core energy Ec and the KTHNY transition: 

Next, we analyze our results in the light of the KTHNY- theory The unrenormalized 
K = IGvr at pc = 0.904 (P^ = 8.92) (see Fig. (2), lower arrow) which implies that a weak first 
order transition from solid to liquid preempts a KTHNY- solid- hexatic transition. However, 
the value of K is renormalized by the presence of dislocations. We can estimate the extent 
of this renormalization from our data. 

The dislocation pair probability 

Prf = exp{-2E,)Z{K) (9) 
where Z{K) is the "internal partition function" of a dislocation pair and is given by [^], 

Where we have set the core radius = a, the lattice parameter. The core energy of a 
dislocation is a difficult quantity to obtain from a simulation, though it has been attempted 
in the past [p5| , p6[ . In our case, an ansatz, which gives excellent results in the 2D- XY- 



model [T^, and identifies the rejection ratio p a.s p = pd can be used to obtain E^, see 



Fig. (7). Throughout the relevant region is safely above the limit E^ > 2.7 p!7| , p!6 
At the transition the Ec ^ Q which is in good agreement the results of Murray and Van 
Winkle |^ {Ef. ~ 5.6) from experiments on 2-d charge stabilised colloids and of Zahn et. 
al. [^] {Ec ~ 4.) for paramagnetic colloids. 
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Finally, to obtain the melting density we use the unrenormalised K and y — exp{—Ec{K)) 
as inputs to the KTHNY recursion relations (Eqs.(4)) and solve them numerically by a 
standard Euler discretization to obtain Kr, see Fig. (8). The melting density obtained from 
our value for Kr is pc — .916 and Pc — 9.39 (Fig. (2), upper arrow). This means that the 
KTHNY- transition now preceeds the first order transition and the solid transforms to the 
hexatic phase. 

IV. SUMMARY AND CONCLUSION 

We have simulated a dislocation free triangular solid of hard disks using a constrained 
Monte Carlo algorithm. Using a block analysis scheme we calculate the finite size scaled 
elastic constants of this solid. From the number of times the system attempts to violate our 
no- dislocation constraint we can obtain (virtual) dislocation probabilities and hence the 
core energy. The absence of a phase transitions in our system implies that all correlation 
lengths remain finite and the problem of slow equilibration of defect densities is eliminated. 
In effect we obtain highly accurate values of the unrenormalized coupling constant K and the 
defect fugacity y which can be used as inputs to the KTHNY recursion relations. Numerical 
solution of these recursion relations then yields the renormalized coupling Kr and hence the 
density and pressure of the solid to hexatic melting transition. 

We can draw a few very precise conclusions from our results. Firstly, a solid without 
dislocations is stable against fluctuations of the amplitude of the solid order parameter 
and against long wavelength phonons. So any melting transition mediated by phonon or 
amplitude fluctuation is ruled out in our system. Secondly, the core energy Ec> 2.7 at the 
transition so KTHNY perturbation theory is valid though numerical values of nonuniversal 
quantities may depend on the order of the perturbation analysis. Thirdly, solution of the 
recursion relations shows that a KTHNY transition at Pc — 9.39 preempts the first order 
transition at Pi — 9.2. Since these transitions, as well as the hexatic -hquid KTHNY 
transition lies so close to each other, the effect of, as yet unknown, higher order corrections 
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to the recursion relations may need to be examined in the future [0. Due to this caveat, 
our conclusion that a hexatic phase exists over some region of density exceeding p = .899 
still must be taken as preliminary. Also, in actual simulations, cross over effects near the 
bicritical point, where two critical lines corresponding to the liquid -hexatic and hexatic -solid 
transitions meet a first order liquid -solid line (see for e.g. Ref. |]2D| for a lattice model where 
such a situation is discussed) may complicate the analysis of the data, which may, in part, 
explain the confusion which persists in the literature on this subject. In systems with softer 
potentials [Q, the signature of a KTHNY transition appears to be more pronounced []3l| . 
In future, we would like to analyze more complicated systems eg. laser induced reentrant 
melting of charge stabilized colloids , and the influence of other defect variables eg. grain 
boundaries p3l on elastic constants and melting behaviour. 
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p 




P 


B 




p X 10^ 


K/IQti 


0.88 


10^ 


8.117 


27.69 


11.63 


0.36823 


0.8550 


0.9 


10^ 


8.777 


32.47 


13.87 


0.20358 


0.9925 


0.905 


10^ 


8.957 


33.67 


14.46 


0.17386 


1.0271 


0.910 


10^ 


9.145 


35.38 


15.22 


0.14469 


1.0744 


0.915 


10^ 


9.342 


37.09 


15.99 


0.11706 


1.1225 


0.920 


10^ 


9.545 


38.48 


16.88 


0.09532 


1.1722 


0.925 


10^ 


9.759 


40.67 


17.88 


0.07513 


1.2337 


0.930 


10^ 


9.982 


42.72 


18.90 


0.05967 


1.2948 


0.935 


10^ 


10.217 


44.69 


19.91 


0.04643 


1.3538 


0.94 


2x10^ 


10.462 


46.85 


21.45 


0.03432 


1.4382 


0.95 


10^ 


10.996 


52.14 


24.10 


0.01855 


1.5945 


0.96 


2x10^ 


11.586 


59.67 


27.61 


0.00901 


1.8067 


0.97 


10^ 


12.251 


67.45 


31.59 


0.00370 


2.0379 


0.98 


2x10^ 


13.003 


79.20 


36.62 


0.00137 


2.3479 


0.99 


10^ 


13.862 


89.98 


42.60 


0.00041 


2.6835 


1. 


49400 


14.843 


104.78 


50.25 


0.00009 


3.1206 


1.02 


10^ 


17.301 


148.88 


69.91 


0.0 


4.2854 


1.04 


10^ 


20.714 


212.02 


102.02 


0.0 


6.0857 


1.06 


10^ 




319.07 


158.69 


0.0 


9.1874 


1.08 


10^ 




531.24 


268.02 


0.0 


15.1567 


1.1 


10^ 




1018.49 


526.94 


0.0 


29.0094 



Table I Pressure P, bulk modulus B, effective sfiear modulus Heff, ratio of moves rejected 
due to the zero dislocation density constraint p and the (unrenormalized) couphng constant 
K/IQt: as a function of the density p. The total number of configurations used for the 
averages Nc is also listed. The pressure P was obtained by integrating B below p — 1.049. 
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FIGURES 



I - ordered 
II - disordered 




1/K 1/16 71 

FIG. 1. Schematic flows of the coupling constant K and the defect fugacity y under the action 
of the KTHNY recursion relations. The dashed line is the separatrix whose intersection with the 
line of initial state (solid line connecting filled circles, y{T,l = 0), K~^(T,l = 0)) determines the 
transition point Tc 
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p 



0.2 0.4 0.6 



1 1.2 




0.86 0.87 0.88 0.89 0.9 0.91 0.92 

P P 
FIG. 2. Equation of state of the hard disk -hquid and -sohd. (a) Liquid: light sohd hne; 



semi-empirical form of Santos et al. |23], O; data from Hoover et al. ||2^. Solid: bold solid line; our 
results, +; data of Wojciechowski and Brahka dotted line; position of coexistence pressure as 
seen in all studies observing a first order phase transition, (b) Expanded view of (a) near the phase 
transition region. O; results of Jaster Q for 128 x 128 particles, +; same for 256 x 256 particles. 
Light solid line; polynomial fit to Jaster 's data, bold solid line; our data for the solid (as in (a)). 
Horizontal dotted line: as in (a), dotted curve; semi-empirical form for the equation of state of the 



hard disk liquid of Santos et. al. |23]. Arrows: lower arrow; position of the KTHNY-transition 
with bare values for K, upper arrow; same with renormalized calculated from our simulations. 
Note that the accuracy of Jaster's data is smaller than the size of the symbols for p < .9, while 
for p > .9 there may be systematic finite size effects and finite observation time effects possibly 
invalidating the data. 
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0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 

Pi 

FIG. 3. The number of hard disks with five fold (ns, <> and light solid line) and seven fold 
(ny, + and bold solid line) co-ordination after 10^ Monte Carlo steps per particle for a iV = 3120 
particle system, plotted against pi (see text). Note that ns ^ wj for pi larger than .92. 
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0.85 0.9 0.95 1 1.05 1.1 1.15 0.85 0.9 0.95 1 1.05 1.1 1.15 

P P 
FIG. 5. Elastic moduli in the thermodynamic limit: (a) bulk B and (b) shear ^eff- ^ - our 

work (error bars are much smaller than the symbol size); + Wojciechowski and Brahka [10|, solid 



line; free volume theory |22], dashed line; polynomial fit given by P. Bladon and D. Frenkel p 
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12.72 



9.15 



9.2 
P 



9.25 



FIG. 6. Gibbs Free energies g{P) as a function of pressure: dotted line; metastable liquid 
using semi empirical form of Santos et. al. |23|; bold solid line; using Jaster's results series of 
light solid lines; gibbs free energy of the solid where we reduced the reference free energy from the 
value quoted by Bladon and Frenkel ||2^ (see text) by 3.3, 3.35 and 3.4%. 
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3 4 
i^/167r 

FIG. 7. Calculated core energy Ec (O) as a function of K/IGtt. The straight line is a linear 



least square fit. Note that Ec> 2.7 throughout. 
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0.905 0.91 0.915 0.92 

P 

FIG. 8. Renormalization of K/lQir vs density p for the hard disk solid. The renormahzed 
Kr/Wtt (bold solid line)i is obtained from the recursion relations Eq. (4) which were solved by 
Euler discretization using a step size SI = 0.001 upto a final I = 100 starting from the initial input 
(light solid line). Dotted line: K = IQ-k. 
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